/* -*- c -*- */

/*
 *****************************************************************************
 **                            INCLUDES                                     **
 *****************************************************************************
 */
#define PY_SSIZE_T_CLEAN
#include <Python.h>

#define NPY_NO_DEPRECATED_API NPY_API_VERSION
#if defined(NPY_INTERNAL_BUILD)
#undef NPY_INTERNAL_BUILD
#endif
// for add_INT32_negative_indexed
#define NPY_TARGET_VERSION NPY_2_1_API_VERSION
#include "numpy/arrayobject.h"
#include "numpy/ufuncobject.h"
#include "numpy/ndarrayobject.h"
#include "numpy/npy_math.h"



#include "npy_config.h"
#include "npy_cpu_features.h"
#include "npy_cpu_dispatch.h"
#include "numpy/npy_cpu.h"
#include "npy_import.h"
#include "numpy/dtype_api.h"

/*
 *****************************************************************************
 **                            BASICS                                       **
 *****************************************************************************
 */

#define INIT_OUTER_LOOP_1       \
    npy_intp dN = *dimensions++;\
    npy_intp N_;                \
    npy_intp s0 = *steps++;

#define INIT_OUTER_LOOP_2       \
    INIT_OUTER_LOOP_1           \
    npy_intp s1 = *steps++;

#define INIT_OUTER_LOOP_3       \
    INIT_OUTER_LOOP_2           \
    npy_intp s2 = *steps++;

#define INIT_OUTER_LOOP_4       \
    INIT_OUTER_LOOP_3           \
    npy_intp s3 = *steps++;

#define BEGIN_OUTER_LOOP_2      \
    for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1) {

#define BEGIN_OUTER_LOOP_3      \
    for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) {

#define BEGIN_OUTER_LOOP_4      \
    for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2, args[3] += s3) {

#define END_OUTER_LOOP  }


/*
 *****************************************************************************
 **                             UFUNC LOOPS                                 **
 *****************************************************************************
 */

static void
always_error_loop(
        char **NPY_UNUSED(args), npy_intp const *NPY_UNUSED(dimensions),
        npy_intp const *NPY_UNUSED(steps), void *NPY_UNUSED(func))
{
    NPY_ALLOW_C_API_DEF
    NPY_ALLOW_C_API;
    PyErr_SetString(PyExc_RuntimeError, "How unexpected :)!");
    NPY_DISABLE_C_API;
    return;
}


char *inner1d_signature = "(i),(i)->()";

/**begin repeat

   #TYPE=INTP,DOUBLE#
   #typ=npy_intp,npy_double#
*/

/*
 *  This implements the function
 *        out[n] = sum_i { in1[n, i] * in2[n, i] }.
 */

static void
@TYPE@_inner1d(char **args, npy_intp const *dimensions, npy_intp const *steps, void *NPY_UNUSED(func))
{
    INIT_OUTER_LOOP_3
    npy_intp di = dimensions[0];
    npy_intp i;
    npy_intp is1=steps[0], is2=steps[1];
    BEGIN_OUTER_LOOP_3
        char *ip1=args[0], *ip2=args[1], *op=args[2];
        @typ@ sum = 0;
        for (i = 0; i < di; i++) {
            sum += (*(@typ@ *)ip1) * (*(@typ@ *)ip2);
            ip1 += is1;
            ip2 += is2;
        }
        *(@typ@ *)op = sum;
    END_OUTER_LOOP
}

/**end repeat**/

char *innerwt_signature = "(i),(i),(i)->()";

/**begin repeat

   #TYPE=INTP,DOUBLE#
   #typ=npy_intp,npy_double#
*/


/*
 *  This implements the function
 *        out[n] = sum_i { in1[n, i] * in2[n, i] * in3[n, i] }.
 */

static void
@TYPE@_innerwt(char **args, npy_intp const *dimensions, npy_intp const *steps, void *NPY_UNUSED(func))
{
    INIT_OUTER_LOOP_4
    npy_intp di = dimensions[0];
    npy_intp i;
    npy_intp is1=steps[0], is2=steps[1], is3=steps[2];
    BEGIN_OUTER_LOOP_4
        char *ip1=args[0], *ip2=args[1], *ip3=args[2], *op=args[3];
        @typ@ sum = 0;
        for (i = 0; i < di; i++) {
            sum += (*(@typ@ *)ip1) * (*(@typ@ *)ip2) * (*(@typ@ *)ip3);
            ip1 += is1;
            ip2 += is2;
            ip3 += is3;
        }
        *(@typ@ *)op = sum;
    END_OUTER_LOOP
}

/**end repeat**/

char *matrix_multiply_signature = "(m,n),(n,p)->(m,p)";
/* for use with matrix_multiply code, but different signature */
char *matmul_signature = "(m?,n),(n,p?)->(m?,p?)";

/**begin repeat

   #TYPE=FLOAT,DOUBLE,INTP#
   #typ=npy_float,npy_double,npy_intp#
*/

/*
 *  This implements the function
 *        out[k, m, p] = sum_n { in1[k, m, n] * in2[k, n, p] }.
 */

static void
@TYPE@_matrix_multiply(char **args, npy_intp const *dimensions, npy_intp const *steps, void *NPY_UNUSED(func))
{
    /* no BLAS is available */
    INIT_OUTER_LOOP_3
    npy_intp dm = dimensions[0];
    npy_intp dn = dimensions[1];
    npy_intp dp = dimensions[2];
    npy_intp m,n,p;
    npy_intp is1_m=steps[0], is1_n=steps[1], is2_n=steps[2], is2_p=steps[3],
         os_m=steps[4], os_p=steps[5];
    npy_intp ib1_n = is1_n*dn;
    npy_intp ib2_n = is2_n*dn;
    npy_intp ib2_p = is2_p*dp;
    npy_intp ob_p  = os_p *dp;
    if (dn == 0) {
        /* No operand, need to zero the output */
        BEGIN_OUTER_LOOP_3
            char *op=args[2];
            for (m = 0; m < dm; m++) {
                for (p = 0; p < dp; p++) {
                    *(@typ@ *)op = 0;
                    op  +=  os_p;
                }
                op  +=  os_m - ob_p;
            }
        END_OUTER_LOOP
        return;
    }
    BEGIN_OUTER_LOOP_3
        char *ip1=args[0], *ip2=args[1], *op=args[2];
        for (m = 0; m < dm; m++) {
            for (n = 0; n < dn; n++) {
                @typ@ val1 = (*(@typ@ *)ip1);
                for (p = 0; p < dp; p++) {
                    if (n == 0) *(@typ@ *)op = 0;
                    *(@typ@ *)op += val1 * (*(@typ@ *)ip2);
                    ip2 += is2_p;
                    op  +=  os_p;
                }
                ip2 -= ib2_p;
                op  -=  ob_p;
                ip1 += is1_n;
                ip2 += is2_n;
            }
            ip1 -= ib1_n;
            ip2 -= ib2_n;
            ip1 += is1_m;
            op  +=  os_m;
        }
    END_OUTER_LOOP
}

/**end repeat**/

char *cross1d_signature = "(3),(3)->(3)";

/**begin repeat

   #TYPE=INTP,DOUBLE#
   #typ=npy_intp, npy_double#
*/

/*
 *  This implements the cross product:
 *        out[n, 0] = in1[n, 1]*in2[n, 2] - in1[n, 2]*in2[n, 1]
 *        out[n, 1] = in1[n, 2]*in2[n, 0] - in1[n, 0]*in2[n, 2]
 *        out[n, 2] = in1[n, 0]*in2[n, 1] - in1[n, 1]*in2[n, 0]
 */
static void
@TYPE@_cross1d(char **args, npy_intp const *dimensions, npy_intp const *steps, void *NPY_UNUSED(func))
{
    INIT_OUTER_LOOP_3
    npy_intp is1=steps[0], is2=steps[1], os = steps[2];
    BEGIN_OUTER_LOOP_3
        @typ@ i1_x = *(@typ@ *)(args[0] + 0*is1);
        @typ@ i1_y = *(@typ@ *)(args[0] + 1*is1);
        @typ@ i1_z = *(@typ@ *)(args[0] + 2*is1);

        @typ@ i2_x = *(@typ@ *)(args[1] + 0*is2);
        @typ@ i2_y = *(@typ@ *)(args[1] + 1*is2);
        @typ@ i2_z = *(@typ@ *)(args[1] + 2*is2);
        char *op = args[2];

        *(@typ@ *)op = i1_y * i2_z - i1_z * i2_y;
        op += os;
        *(@typ@ *)op = i1_z * i2_x - i1_x * i2_z;
        op += os;
        *(@typ@ *)op = i1_x * i2_y - i1_y * i2_x;
    END_OUTER_LOOP
}

/**end repeat**/

char *euclidean_pdist_signature = "(n,d)->(p)";

/**begin repeat

   #TYPE=FLOAT,DOUBLE#
   #typ=npy_float,npy_double#
   #sqrt_func=sqrtf,sqrt#
*/

/*
 *  This implements the function
 *        out[j*(2*n-3-j)+k-1] = sum_d { (in1[j, d] - in1[k, d])^2 }
 *  with 0 < k < j < n, i.e. computes all unique pairwise euclidean distances.
 */

static void
@TYPE@_euclidean_pdist(char **args, npy_intp const *dimensions, npy_intp const *steps,
                       void *NPY_UNUSED(func))
{
    INIT_OUTER_LOOP_2
    npy_intp len_n = *dimensions++;
    npy_intp len_d = *dimensions++;
    npy_intp stride_n = *steps++;
    npy_intp stride_d = *steps++;
    npy_intp stride_p = *steps;

    assert(len_n * (len_n - 1) / 2 == *dimensions);

    BEGIN_OUTER_LOOP_2
        const char *data_this = (const char *)args[0];
        char *data_out = args[1];
        npy_intp n;
        for (n = 0; n < len_n; ++n) {
            const char *data_that = data_this + stride_n;
            npy_intp nn;
            for (nn = n + 1; nn < len_n; ++nn) {
                const char *ptr_this = data_this;
                const char *ptr_that = data_that;
                @typ@ out = 0;
                npy_intp d;
                for (d = 0; d < len_d; ++d) {
                    const @typ@ delta = *(const @typ@ *)ptr_this -
                                        *(const @typ@ *)ptr_that;
                    out += delta * delta;
                    ptr_this += stride_d;
                    ptr_that += stride_d;
                }
                *(@typ@ *)data_out = @sqrt_func@(out);
                data_that += stride_n;
                data_out += stride_p;
            }
            data_this += stride_n;
        }
    END_OUTER_LOOP
}

/**end repeat**/

char *cumsum_signature = "(i)->(i)";

/*
 *  This implements the function
 *        out[n] = sum_i^n in[i]
 */

/**begin repeat

   #TYPE=INTP,DOUBLE#
   #typ=npy_intp,npy_double#
*/

static void
@TYPE@_cumsum(char **args, npy_intp const *dimensions, npy_intp const *steps, void *NPY_UNUSED(func))
{
    INIT_OUTER_LOOP_2
    npy_intp di = dimensions[0];
    npy_intp i;
    npy_intp is=steps[0], os=steps[1];
    BEGIN_OUTER_LOOP_2
        char *ip=args[0], *op=args[1];
        @typ@ cumsum = 0;
        for (i = 0; i < di; i++, ip += is, op += os) {
            cumsum += (*(@typ@ *)ip);
            *(@typ@ *)op = cumsum;
        }
    END_OUTER_LOOP
}

/**end repeat**/

static int
INT32_negative(PyArrayMethod_Context *NPY_UNUSED(context),
               char **args, npy_intp const *dimensions,
               npy_intp const *steps, NpyAuxData *NPY_UNUSED(func))
{
    npy_intp di = dimensions[0];
    npy_intp i;
    npy_intp is=steps[0], os=steps[1];
    char *ip=args[0], *op=args[1];
    for (i = 0; i < di; i++, ip += is, op += os) {
        if (i == 3) {
            *(int32_t *)op = - 100;
        } else {
            *(int32_t *)op = - *(int32_t *)ip;
        }
    }
    return 0;
}


static int
INT32_negative_indexed(PyArrayMethod_Context *NPY_UNUSED(context),
                           char * const*args, npy_intp const *dimensions,
                           npy_intp const *steps, NpyAuxData *NPY_UNUSED(func))
{
    char *ip1 = args[0];
    char *indxp = args[1];
    npy_intp is1 = steps[0], isindex = steps[1];
    npy_intp n = dimensions[0];
    npy_intp shape = steps[3];
    npy_intp i;
    int32_t *indexed;
    for(i = 0; i < n; i++, indxp += isindex) {
        npy_intp indx = *(npy_intp *)indxp;
        if (indx < 0) {
            indx += shape;
        }
        indexed = (int32_t *)(ip1 + is1 * indx);
        if (i == 3) {
            *indexed = -200;
        } else {
            *indexed = - *indexed;
        }
    }
    return 0;
}



/*  The following lines were generated using a slightly modified
    version of code_generators/generate_umath.py and adding these
    lines to defdict:

defdict = {
'inner1d' :
    Ufunc(2, 1, None_,
        r'''inner on the last dimension and broadcast on the rest \n"
        "     \"(i),(i)->()\" \n''',
        TD('ld'),
        ),
'innerwt' :
    Ufunc(3, 1, None_,
        r'''inner1d with a weight argument \n"
        "     \"(i),(i),(i)->()\" \n''',
        TD('ld'),
        ),
}

*/

static PyUFuncGenericFunction always_error_functions[] = { always_error_loop };
static void *const always_error_data[] = { (void *)NULL };
static const char always_error_signatures[] = { NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE };
static PyUFuncGenericFunction inner1d_functions[] = { INTP_inner1d, DOUBLE_inner1d };
static void *const inner1d_data[] = { (void *)NULL, (void *)NULL };
static const char inner1d_signatures[] = { NPY_INTP, NPY_INTP, NPY_INTP, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE };
static PyUFuncGenericFunction innerwt_functions[] = { INTP_innerwt, DOUBLE_innerwt };
static void *const innerwt_data[] = { (void *)NULL, (void *)NULL };
static const char innerwt_signatures[] = { NPY_INTP, NPY_INTP, NPY_INTP, NPY_INTP, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE };
static PyUFuncGenericFunction matrix_multiply_functions[] = { INTP_matrix_multiply, FLOAT_matrix_multiply, DOUBLE_matrix_multiply };
static void *const matrix_multiply_data[] = { (void *)NULL, (void *)NULL, (void *)NULL };
static const char matrix_multiply_signatures[] = { NPY_INTP, NPY_INTP, NPY_INTP,  NPY_FLOAT, NPY_FLOAT, NPY_FLOAT,  NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE };
static PyUFuncGenericFunction cross1d_functions[] = { INTP_cross1d, DOUBLE_cross1d };
static void *const cross1d_data[] = { (void *)NULL, (void *)NULL };
static const char cross1d_signatures[] = { NPY_INTP, NPY_INTP, NPY_INTP, NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE };
static PyUFuncGenericFunction euclidean_pdist_functions[] =
                            { FLOAT_euclidean_pdist, DOUBLE_euclidean_pdist };
static void *const eucldiean_pdist_data[] = { (void *)NULL, (void *)NULL };
static const char euclidean_pdist_signatures[] = { NPY_FLOAT, NPY_FLOAT,
                                             NPY_DOUBLE, NPY_DOUBLE };

static PyUFuncGenericFunction cumsum_functions[] = { INTP_cumsum, DOUBLE_cumsum };
static void *const cumsum_data[] = { (void *)NULL, (void *)NULL };
static const char cumsum_signatures[] = { NPY_INTP, NPY_INTP, NPY_DOUBLE, NPY_DOUBLE };


static int
addUfuncs(PyObject *dictionary) {
    PyObject *f;

    f = PyUFunc_FromFuncAndData(always_error_functions, always_error_data,
            always_error_signatures, 1, 2, 1, PyUFunc_None, "always_error",
            "simply, broken, ufunc that sets an error (but releases the GIL).",
            0);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "always_error", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(always_error_functions,
            always_error_data, always_error_signatures, 1, 2, 1, PyUFunc_None,
            "always_error_gufunc",
            "simply, broken, gufunc that sets an error (but releases the GIL).",
            0, "(i),()->()");
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "always_error_gufunc", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(inner1d_functions, inner1d_data,
                    inner1d_signatures, 2, 2, 1, PyUFunc_None, "inner1d",
                    "inner on the last dimension and broadcast on the rest \n"
                    "     \"(i),(i)->()\" \n",
                    0, inner1d_signature);
    /*
     * yes, this should not happen, but I (MHvK) just spent an hour looking at
     * segfaults because I screwed up something that seemed totally unrelated.
     */
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "inner1d", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(innerwt_functions, innerwt_data,
                    innerwt_signatures, 2, 3, 1, PyUFunc_None, "innerwt",
                    "inner1d with a weight argument \n"
                    "     \"(i),(i),(i)->()\" \n",
                    0, innerwt_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "innerwt", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(matrix_multiply_functions,
                    matrix_multiply_data, matrix_multiply_signatures,
                    3, 2, 1, PyUFunc_None, "matrix_multiply",
                    "matrix multiplication on last two dimensions \n"
                    "     \"(m,n),(n,p)->(m,p)\" \n",
                    0, matrix_multiply_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "matrix_multiply", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(matrix_multiply_functions,
                    matrix_multiply_data, matrix_multiply_signatures,
                    3, 2, 1, PyUFunc_None, "matmul",
                    "matmul on last two dimensions, with some being optional\n"
                    "     \"(m?,n),(n,p?)->(m?,p?)\" \n",
                    0, matmul_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "matmul", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(euclidean_pdist_functions,
                    eucldiean_pdist_data, euclidean_pdist_signatures,
                    2, 1, 1, PyUFunc_None, "euclidean_pdist",
                    "pairwise euclidean distance on last two dimensions \n"
                    "     \"(n,d)->(p)\" \n",
                    0, euclidean_pdist_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "euclidean_pdist", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(cumsum_functions,
                    cumsum_data, cumsum_signatures,
                    2, 1, 1, PyUFunc_None, "cumsum",
                    "Cumulative sum of the input (n)->(n)\n",
                    0, cumsum_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "cumsum", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(inner1d_functions, inner1d_data,
                    inner1d_signatures, 2, 2, 1, PyUFunc_None, "inner1d_no_doc",
                    NULL,
                    0, inner1d_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "inner1d_no_doc", f);
    Py_DECREF(f);
    f = PyUFunc_FromFuncAndDataAndSignature(cross1d_functions, cross1d_data,
                    cross1d_signatures, 2, 2, 1, PyUFunc_None, "cross1d",
                    "cross product on the last dimension and broadcast on the rest \n"\
                    "     \"(3),(3)->(3)\" \n",
                    0, cross1d_signature);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "cross1d", f);
    Py_DECREF(f);

    f = PyUFunc_FromFuncAndDataAndSignature(NULL, NULL,
            NULL, 0, 0, 0, PyUFunc_None, "_pickleable_module_global.ufunc",
            "A dotted name for pickle testing, does nothing.", 0, NULL);
    if (f == NULL) {
        return -1;
    }
    PyDict_SetItemString(dictionary, "_pickleable_module_global_ufunc", f);
    Py_DECREF(f);

    return 0;
}


static PyObject *
UMath_Tests_test_signature(PyObject *NPY_UNUSED(dummy), PyObject *args)
{
    int nin, nout, i;
    PyObject *signature=NULL, *sig_str=NULL;
    PyUFuncObject *f=NULL;
    PyObject *core_num_dims=NULL, *core_dim_ixs=NULL;
    PyObject *core_dim_flags=NULL, *core_dim_sizes=NULL;
    int core_enabled;
    int core_num_ixs = 0;

    if (!PyArg_ParseTuple(args, "iiO", &nin, &nout, &signature)) {
        return NULL;
    }

    if (PyBytes_Check(signature)) {
        sig_str = signature;
    } else if (PyUnicode_Check(signature)) {
        sig_str = PyUnicode_AsUTF8String(signature);
    } else {
        PyErr_SetString(PyExc_ValueError, "signature should be a string");
        return NULL;
    }

    f = (PyUFuncObject*)PyUFunc_FromFuncAndDataAndSignature(
        NULL, NULL, NULL,
        0, nin, nout, PyUFunc_None, "no name",
        "doc:none",
        1, PyBytes_AS_STRING(sig_str));
    if (sig_str != signature) {
        Py_DECREF(sig_str);
    }
    if (f == NULL) {
        return NULL;
    }
    core_enabled = f->core_enabled;
    /*
     * Don't presume core_num_dims and core_dim_ixs are defined;
     * they currently are even if core_enabled=0, but there's no real
     * reason they should be.  So avoid segfaults if we change our mind.
     */
    if (f->core_num_dims != NULL) {
        core_num_dims = PyTuple_New(f->nargs);
        if (core_num_dims == NULL) {
            goto fail;
        }
        for (i = 0; i < f->nargs; i++) {
            PyObject *val = PyLong_FromLong(f->core_num_dims[i]);
            PyTuple_SET_ITEM(core_num_dims, i, val);
            core_num_ixs += f->core_num_dims[i];
        }
    }
    else {
        Py_INCREF(Py_None);
        core_num_dims = Py_None;
    }
    if (f->core_dim_ixs != NULL) {
        core_dim_ixs = PyTuple_New(core_num_ixs);
        if (core_dim_ixs == NULL) {
            goto fail;
        }
        for (i = 0; i < core_num_ixs; i++) {
            PyObject *val = PyLong_FromLong(f->core_dim_ixs[i]);
            PyTuple_SET_ITEM(core_dim_ixs, i, val);
        }
    }
    else {
        Py_INCREF(Py_None);
        core_dim_ixs = Py_None;
    }
    if (f->core_dim_flags != NULL) {
        core_dim_flags = PyTuple_New(f->core_num_dim_ix);
        if (core_dim_flags == NULL) {
            goto fail;
        }
        for (i = 0; i < f->core_num_dim_ix; i++) {
            PyObject *val = PyLong_FromLong(f->core_dim_flags[i]);
            PyTuple_SET_ITEM(core_dim_flags, i, val);
        }
    }
    else {
        Py_INCREF(Py_None);
        core_dim_flags = Py_None;
    }
    if (f->core_dim_sizes != NULL) {
        core_dim_sizes = PyTuple_New(f->core_num_dim_ix);
        if (core_dim_sizes == NULL) {
            goto fail;
        }
        for (i = 0; i < f->core_num_dim_ix; i++) {
            PyObject *val = PyLong_FromLong(f->core_dim_sizes[i]);
            PyTuple_SET_ITEM(core_dim_sizes, i, val);
        }
    }
    else {
        Py_INCREF(Py_None);
        core_dim_sizes = Py_None;
    }
    Py_DECREF(f);
    return Py_BuildValue("iNNNN", core_enabled, core_num_dims,
                         core_dim_ixs, core_dim_flags, core_dim_sizes);

fail:
    Py_XDECREF(f);
    Py_XDECREF(core_num_dims);
    Py_XDECREF(core_dim_ixs);
    Py_XDECREF(core_dim_flags);
    Py_XDECREF(core_dim_sizes);
    return NULL;
}

// Testing the utilities of the CPU dispatcher
#include "_umath_tests.dispatch.h"
NPY_CPU_DISPATCH_DECLARE(extern const char *_umath_tests_dispatch_var)
NPY_CPU_DISPATCH_DECLARE(const char *_umath_tests_dispatch_func, (void))
NPY_CPU_DISPATCH_DECLARE(void _umath_tests_dispatch_attach, (PyObject *list))

static PyObject *
UMath_Tests_test_dispatch(PyObject *NPY_UNUSED(dummy), PyObject *NPY_UNUSED(dummy2))
{
    const char *highest_func, *highest_var;
    NPY_CPU_DISPATCH_CALL(highest_func = _umath_tests_dispatch_func, ());
    NPY_CPU_DISPATCH_CALL(highest_var  = _umath_tests_dispatch_var);
    const char *highest_func_xb = "nobase", *highest_var_xb = "nobase";
    NPY_CPU_DISPATCH_CALL_XB(highest_func_xb = _umath_tests_dispatch_func, ());
    NPY_CPU_DISPATCH_CALL_XB(highest_var_xb  = _umath_tests_dispatch_var);

    PyObject *dict = PyDict_New(), *item;
    if (dict == NULL) {
        return NULL;
    }
    /**begin repeat
     * #str = func, var, func_xb, var_xb#
    */
    item = PyUnicode_FromString(highest_@str@);
    if (item == NULL || PyDict_SetItemString(dict, "@str@", item) < 0) {
        goto err;
    }
    Py_DECREF(item);
    /**end repeat**/
    item = PyList_New(0);
    if (item == NULL || PyDict_SetItemString(dict, "all", item) < 0) {
        goto err;
    }
    NPY_CPU_DISPATCH_CALL_ALL(_umath_tests_dispatch_attach, (item));
    Py_SETREF(item, NULL);
    if (PyErr_Occurred()) {
        goto err;
    }
    return dict;
err:
    Py_XDECREF(item);
    Py_DECREF(dict);
    return NULL;
}

static int
add_INT32_negative_indexed(PyObject *module, PyObject *dict) {
    PyObject * negative = PyUFunc_FromFuncAndData(NULL, NULL, NULL, 0, 1, 1,
                                    PyUFunc_Zero, "indexed_negative", NULL, 0);
    if (negative == NULL) {
        return -1;
    }
    PyArray_DTypeMeta *dtypes[] = {&PyArray_Int32DType, &PyArray_Int32DType};

    PyType_Slot slots[] = {
        {NPY_METH_contiguous_indexed_loop, INT32_negative_indexed},
        {NPY_METH_strided_loop, INT32_negative},
        {0, NULL}
    };

    PyArrayMethod_Spec spec = {
        .name = "negative_indexed_loop",
        .nin = 1,
        .nout = 1,
        .dtypes = dtypes,
        .slots = slots,
        .flags = NPY_METH_NO_FLOATINGPOINT_ERRORS
    };

    if (PyUFunc_AddLoopFromSpec(negative, &spec) < 0) {
        Py_DECREF(negative);
        return -1;
    }
    PyDict_SetItemString(dict, "indexed_negative", negative);
    Py_DECREF(negative);
    return 0;
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
// Define the gufunc 'conv1d_full'
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define MAX(a, b) (((a) < (b)) ? (b) : (a))

int conv1d_full_process_core_dims(PyUFuncObject *ufunc,
                                  npy_intp *core_dim_sizes)
{
    //
    // core_dim_sizes will hold the core dimensions [m, n, p].
    // p will be -1 if the caller did not provide the out argument.
    //
    npy_intp m = core_dim_sizes[0];
    npy_intp n = core_dim_sizes[1];
    npy_intp p = core_dim_sizes[2];
    npy_intp required_p = m + n - 1;

    if (m == 0 && n == 0) {
        PyErr_SetString(PyExc_ValueError,
            "conv1d_full: both inputs have core dimension 0; the function "
            "requires that at least one input has positive size.");
        return -1;
    }
    if (p == -1) {
        core_dim_sizes[2] = required_p;
        return 0;
    }
    if (p != required_p) {
        PyErr_Format(PyExc_ValueError,
                "conv1d_full: the core dimension p of the out parameter "
                "does not equal m + n - 1, where m and n are the core "
                "dimensions of the inputs x and y; got m=%zd and n=%zd so "
                "p must be %zd, but got p=%zd.",
                m, n, required_p, p);
        return -1;
    }
    return 0;
}

static void
conv1d_full_double_loop(char **args,
                        npy_intp const *dimensions,
                        npy_intp const *steps,
                        void *NPY_UNUSED(func))
{
    // Input and output arrays
    char *p_x = args[0];
    char *p_y = args[1];
    char *p_out = args[2];
    // Number of loops of pdist calculations to execute.
    npy_intp nloops = dimensions[0];
    // Core dimensions
    npy_intp m = dimensions[1];
    npy_intp n = dimensions[2];
    npy_intp p = dimensions[3]; // Must be m + n - 1.
    // Core strides
    npy_intp x_stride = steps[0];
    npy_intp y_stride = steps[1];
    npy_intp out_stride = steps[2];
    // Inner strides
    npy_intp x_inner_stride = steps[3];
    npy_intp y_inner_stride = steps[4];
    npy_intp out_inner_stride = steps[5];

    for (npy_intp loop = 0; loop < nloops; ++loop, p_x += x_stride,
                                                   p_y += y_stride,
                                                   p_out += out_stride) {
        // Basic implementation of 1d convolution
        for (npy_intp k = 0; k < p; ++k) {
            double sum = 0.0;
            for (npy_intp i = MAX(0, k - n + 1); i < MIN(m, k + 1); ++i) {
                double x_i = *(double *)(p_x + i*x_inner_stride);
                double y_k_minus_i = *(double *)(p_y + (k - i)*y_inner_stride);
                sum +=  x_i * y_k_minus_i;
            }
            *(double *)(p_out + k*out_inner_stride) = sum;
        }
    }
}

static PyUFuncGenericFunction conv1d_full_functions[] = {
    (PyUFuncGenericFunction) &conv1d_full_double_loop
};
static void *const conv1d_full_data[] = {NULL};
static const char conv1d_full_typecodes[] = {NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE};


static PyMethodDef UMath_TestsMethods[] = {
    {"test_signature",  UMath_Tests_test_signature, METH_VARARGS,
     "Test signature parsing of ufunc. \n"
     "Arguments: nin nout signature \n"
     "If fails, it returns NULL. Otherwise it returns a tuple of ufunc "
     "internals. \n",
     },
    {"test_dispatch", UMath_Tests_test_dispatch, METH_NOARGS, NULL},
    {NULL, NULL, 0, NULL}        /* Sentinel */
};

static struct PyModuleDef moduledef = {
        PyModuleDef_HEAD_INIT,
        "_umath_tests",
        NULL,
        -1,
        UMath_TestsMethods,
        NULL,
        NULL,
        NULL,
        NULL
};

/* Initialization function for the module */
PyMODINIT_FUNC PyInit__umath_tests(void) {
    PyObject *m;
    PyObject *d;
    PyObject *version;

    // Initialize CPU features
    if (npy_cpu_init() < 0) {
        return NULL;
    }

    m = PyModule_Create(&moduledef);
    if (m == NULL) {
        return NULL;
    }

    if (PyArray_ImportNumPyAPI() < 0) {
        return NULL;
    }
    if (PyUFunc_ImportUFuncAPI() < 0) {
        return NULL;
    }

    d = PyModule_GetDict(m);

    version = PyUnicode_FromString("0.1");
    PyDict_SetItemString(d, "__version__", version);
    Py_DECREF(version);

    /* Load the ufunc operators into the module's namespace */
    if (addUfuncs(d) < 0) {
        Py_DECREF(m);
        PyErr_Print();
        PyErr_SetString(PyExc_RuntimeError,
                        "cannot load _umath_tests module.");
        return NULL;
    }

    if (add_INT32_negative_indexed(m, d) < 0) {
        Py_DECREF(m);
        PyErr_Print();
        PyErr_SetString(PyExc_RuntimeError,
                        "cannot load _umath_tests module.");
        return NULL;
    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    // Define the gufunc 'conv1d_full'
    // Shape signature is (m),(n)->(p) where p must be m + n - 1.
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    PyUFuncObject *gufunc = (PyUFuncObject *) PyUFunc_FromFuncAndDataAndSignature(
                                conv1d_full_functions,
                                conv1d_full_data,
                                conv1d_full_typecodes,
                                1, 2, 1, PyUFunc_None, "conv1d_full",
                                "convolution of x and y ('full' mode)",
                                0, "(m),(n)->(p)");
    if (gufunc == NULL) {
        Py_DECREF(m);
        return NULL;
    }
    gufunc->process_core_dims_func = &conv1d_full_process_core_dims;

    int status = PyModule_AddObject(m, "conv1d_full", (PyObject *) gufunc);
    if (status == -1) {
        Py_DECREF(gufunc);
        Py_DECREF(m);
        return NULL;
    }

    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

#ifdef Py_GIL_DISABLED
    // signal this module supports running with the GIL disabled
    PyUnstable_Module_SetGIL(m, Py_MOD_GIL_NOT_USED);
#endif

    return m;
}
